// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtSimilarityTransformFilter.C
// ************************************************************************* //

#include <avtSimilarityTransformFilter.h>

#include <avtExtents.h>
#include <vtkDataArray.h>
#include <vtkPointData.h>
#include <vtkRectilinearGrid.h>

#include <BadVectorException.h>


// ****************************************************************************
//  Method: avtSimilarityTransformFilter constructor
//
//  Programmer: meredith -- generated by xml2info
//  Creation:   Mon Sep 24 13:29:16 PST 2001
//
//  Modifications:
//    Kathleen Bonnell, Thu Apr 10 11:07:48 PDT 2003
//    Initialize invM.
//
// ****************************************************************************

avtSimilarityTransformFilter::avtSimilarityTransformFilter()
{
    M = NULL;
    invM = NULL;
}


// ****************************************************************************
//  Method: avtSimilarityTransformFilter destructor
//
//  Programmer: meredith -- generated by xml2info
//  Creation:   Mon Sep 24 13:29:16 PST 2001
//
//  Modifications:
//    Kathleen Bonnell, Thu Apr 10 11:07:48 PDT 2003
//    Delete invM.
//
// ****************************************************************************

avtSimilarityTransformFilter::~avtSimilarityTransformFilter()
{
    if (M)
        M->Delete();
    M = NULL;
    if (invM)
        invM->Delete();
    invM = NULL;
}


// ****************************************************************************
//  Method:      avtSimilarityTransformFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: meredith -- generated by xml2info
//  Creation:   Mon Sep 24 13:29:16 PST 2001
//
//  Modifications:
//    Kathleen Bonnell, Thu Apr 10 11:07:48 PDT 2003
//    Delete invM.
//
//    Kathleen Bonnell, Wed May 21 11:38:23 PDT 2003   
//    Check for bad axis of rotation. 
//
//    Dave Pugmire, Fri May 14 08:04:43 EDT 2010
//    Flag for vector transformations.
//
// ****************************************************************************

void
avtSimilarityTransformFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const SimilarityTransformAttributes*)a;
    if (M)
        M->Delete();
    M = NULL;
    if (invM)
        invM->Delete();
    invM = NULL;

    SetVectorTransform(atts.GetTransformVectors());

    const double *axis = atts.GetRotateAxis();
    if (axis[0] == 0. && axis[1] == 0. && axis[2] == 0.)
    {
        EXCEPTION1(BadVectorException, "Rotation Axis");
        return;
    }
}


// ****************************************************************************
//  Method:  Jeremy Meredith, Mon Sep 24 14:24:14 PDT 2001
//
//  Purpose:
//    Setup the vtk matrix from the transform attributs
//
//  Programmer:  Jeremy Meredith
//  Creation:    September 24, 2001
//
//  Modifications:
//    Kathleen Bonnell, Thu Apr 10 11:07:48 PDT 2003
//    Compute the inverse matrix. 
//
//    Kathleen Biagas, Mon Aug 15 14:17:39 PDT 2016
//    VTK-8, vtkMatrix4x4 no longer has operator[], use SetElement instead.
//
// ****************************************************************************

void
avtSimilarityTransformFilter::SetupMatrix()
{
    if (M)
        return;

    M = vtkMatrix4x4::New();
    M->Identity();
    invM = vtkMatrix4x4::New();
    invM->Identity();

    //
    //  A place to store the inverse transforms for each type.
    //
    vtkMatrix4x4 *IR = vtkMatrix4x4::New(); // inverse Rotate
    IR->Identity();
    vtkMatrix4x4 *IS = vtkMatrix4x4::New(); // inverse Scale
    IS->Identity();
    vtkMatrix4x4 *IT = vtkMatrix4x4::New(); // inverse Transform
    IT->Identity();

    //
    // Do the rotation
    //
    if (atts.GetDoRotate())
    {
        double oX = atts.GetRotateOrigin()[0];
        double oY = atts.GetRotateOrigin()[1];
        double oZ = atts.GetRotateOrigin()[2];

        double X = atts.GetRotateAxis()[0];
        double Y = atts.GetRotateAxis()[1];
        double Z = atts.GetRotateAxis()[2];
        double p = atts.GetRotateAmount();

        // Convert to radians
        if (atts.GetRotateType() == atts.Deg)
            p *= 3.1415926535898/180.;

        // Normalize the axis vector
        double len = sqrt(X*X + Y*Y + Z*Z);
        if (len)
        {
            X /= len;
            Y /= len;
            Z /= len;
        }

        // Compose the quaternion
        double sin_p = sin(p / 2.);
        double cos_p = cos(p / 2.);

        double q[4] = { X*sin_p, Y*sin_p, Z*sin_p, cos_p };

        vtkMatrix4x4 *T1  = vtkMatrix4x4::New();
        vtkMatrix4x4 *R   = vtkMatrix4x4::New();
        vtkMatrix4x4 *T2  = vtkMatrix4x4::New();
        vtkMatrix4x4 *tmp = vtkMatrix4x4::New();

        // Pre-translate
        T1->Identity();
        T1->SetElement(0, 3, -oX);
        T1->SetElement(1, 3, -oY);
        T1->SetElement(2, 3, -oZ);
        
        // Rotate
        R->Identity();
        R->SetElement(0, 0, 1.0 - 2.0 * (q[1] * q[1] + q[2] * q[2]));
        R->SetElement(0, 1, 2.0 * (q[0] * q[1] - q[2] * q[3]));
        R->SetElement(0, 2, 2.0 * (q[2] * q[0] + q[1] * q[3]));

        R->SetElement(1, 0, 2.0 * (q[0] * q[1] + q[2] * q[3]));
        R->SetElement(1, 1, 1.0 - 2.0 * (q[2] * q[2] + q[0] * q[0]));
        R->SetElement(1, 2, 2.0 * (q[1] * q[2] - q[0] * q[3]));

        R->SetElement(2, 0, 2.0 * (q[2] * q[0] - q[1] * q[3]));
        R->SetElement(2, 1, 2.0 * (q[1] * q[2] + q[0] * q[3]));
        R->SetElement(2, 2, 1.0 - 2.0 * (q[1] * q[1] + q[0] * q[0]));

        // Post-translate
        T2->Identity();
        T2->SetElement(0, 3, oX);
        T2->SetElement(1, 3, oY);
        T2->SetElement(2, 3, oZ);

        // apply it
        vtkMatrix4x4::Multiply4x4(T1,  M,  tmp);   M->DeepCopy(tmp);
        vtkMatrix4x4::Multiply4x4(R,   M,  tmp);   M->DeepCopy(tmp);
        vtkMatrix4x4::Multiply4x4(T2,  M,  tmp);   M->DeepCopy(tmp);

        // Create the inverse rotation.
        vtkMatrix4x4 *R_t= vtkMatrix4x4::New();
        vtkMatrix4x4::Transpose(R, R_t);

        vtkMatrix4x4::Multiply4x4(T1,  IR,  tmp);   IR->DeepCopy(tmp);
        vtkMatrix4x4::Multiply4x4(R_t, IR,  tmp);   IR->DeepCopy(tmp);
        vtkMatrix4x4::Multiply4x4(T2,  IR,  tmp);   IR->DeepCopy(tmp);

        R_t->Delete();
        T1 ->Delete();
        R  ->Delete();
        T2 ->Delete();
        tmp->Delete();
    }
    //
    // Do the scale
    //
    if (atts.GetDoScale())
    {
        double oX = atts.GetScaleOrigin()[0];
        double oY = atts.GetScaleOrigin()[1];
        double oZ = atts.GetScaleOrigin()[2];

        double X = atts.GetScaleX();
        double Y = atts.GetScaleY();
        double Z = atts.GetScaleZ();

        vtkMatrix4x4 *T1  = vtkMatrix4x4::New();
        vtkMatrix4x4 *S   = vtkMatrix4x4::New();
        vtkMatrix4x4 *T2  = vtkMatrix4x4::New();
        vtkMatrix4x4 *tmp = vtkMatrix4x4::New();

        // Pre-translate
        T1->Identity();
        T1->SetElement(0, 3, -oX);
        T1->SetElement(1, 3, -oY);
        T1->SetElement(2, 3, -oZ);

        // Scale
        S->Identity();
        S->SetElement(0, 0, X);
        S->SetElement(1, 1, Y);
        S->SetElement(2, 2, Z);

        // Post-translate
        T2->Identity();
        T2->SetElement(0, 3, oX);
        T2->SetElement(1, 3, oY);
        T2->SetElement(2, 3, oZ);

        // apply it
        vtkMatrix4x4::Multiply4x4(T1,  M,  tmp);   M->DeepCopy(tmp);
        vtkMatrix4x4::Multiply4x4(S,   M,  tmp);   M->DeepCopy(tmp);
        vtkMatrix4x4::Multiply4x4(T2,  M,  tmp);   M->DeepCopy(tmp);

        // create the inverse scale
        vtkMatrix4x4 *S_i= vtkMatrix4x4::New();
        S_i->Identity();
        S_i->SetElement(0, 0, (X != 0. ? 1./X : 0.));
        S_i->SetElement(1, 1, (Y != 0. ? 1./Y : 0.));
        S_i->SetElement(2, 2, (Z != 0. ? 1./Z : 0.));
        vtkMatrix4x4::Multiply4x4(T1,  IS,  tmp);   IS->DeepCopy(tmp);
        vtkMatrix4x4::Multiply4x4(S_i, IS,  tmp);   IS->DeepCopy(tmp);
        vtkMatrix4x4::Multiply4x4(T2,  IS,  tmp);   IS->DeepCopy(tmp);

        S_i->Delete(); 
        T1 ->Delete();
        S  ->Delete();
        T2 ->Delete();
        tmp->Delete();
    }
    //
    // Do the translation
    //
    if (atts.GetDoTranslate())
    {
        vtkMatrix4x4 *T = vtkMatrix4x4::New();
        vtkMatrix4x4 *tmp = vtkMatrix4x4::New();

        // Translate
        T->SetElement(0, 3, atts.GetTranslateX());
        T->SetElement(1, 3, atts.GetTranslateY());
        T->SetElement(2, 3, atts.GetTranslateZ());

        // apply it
        vtkMatrix4x4::Multiply4x4(T,  M,   tmp);   M->DeepCopy(tmp);

        // Inverse translation
        IT->SetElement(0, 3, -atts.GetTranslateX());
        IT->SetElement(1, 3, -atts.GetTranslateY());
        IT->SetElement(2, 3, -atts.GetTranslateZ());

        T  ->Delete();
        tmp->Delete();
    }

    //
    // Compose the separate inverse matrices into one.
    //  IT * IS * IR
    // 
    vtkMatrix4x4 *tmp = vtkMatrix4x4::New();
    vtkMatrix4x4::Multiply4x4(IT, invM, tmp); invM->DeepCopy(tmp);
    vtkMatrix4x4::Multiply4x4(IS, invM, tmp); invM->DeepCopy(tmp);
    vtkMatrix4x4::Multiply4x4(IR, invM, tmp); invM->DeepCopy(tmp);

    tmp->Delete();

    IR  ->Delete();
    IS  ->Delete();
    IT  ->Delete();
}


// ****************************************************************************
//  Method: avtSimilarityTransformFilter::PerformRestriciton
//
//  Purpose:
//    Turn on Zone numbers flag if needed, so that original cell array
//    will be propagated throught the pipeline.
//
//  Programmer: Kathleen Bonnell
//  Creation:   November 28, 2001
//
//  Modifications:
//    Kathleen Bonnell, Wed Jun 19 12:28:10 PDT 2002
//    Don't turn off Zone numbers if they have been turned on elsewhere in
//    the pipeline.
//
//    Kathleen Bonnell, Wed Jun 19 13:42:37 PDT 2002
//    Completely removed the code turning off zone numbers.  Why set a flag
//    to false if it is already false?  False is the default setting.
//
//    Kathleen Bonnell, Wed Jun  2 09:08:35 PDT 2004 
//    Added code to turn on node numbers when required. 
//
// ****************************************************************************

avtContract_p
avtSimilarityTransformFilter::ModifyContract(avtContract_p spec)
{
    avtContract_p rv = new avtContract(spec);
    if (rv->GetDataRequest()->MayRequireZones())
    {
        rv->GetDataRequest()->TurnZoneNumbersOn();
    }
    if (rv->GetDataRequest()->MayRequireNodes())
    {
        rv->GetDataRequest()->TurnNodeNumbersOn();
    }
    return rv;
}


// ****************************************************************************
//  Method: avtSimilarityTransformFilter::UpdateDataObjectInfo
//
//  Purpose:
//      If a 2D plot is being revolved into 3D space, then indicate that it is
//      a 3D plot.
//
//  Programmer: Hank Childs
//  Creation:   March 7, 2003
//
//  Modifications:
//  
//    Hank Childs, Thu Jun 17 10:44:20 PDT 2004
//    Set the topological dimension to be higher so that we know we need
//    to apply the facelist filter.  Also transform the extents.
//    
//    Hank Childs, Thu Jul 28 09:07:31 PDT 2005
//    Fix memory leak.
//
//    Hank Childs, Thu Aug 26 13:47:30 PDT 2010
//    Change extents names.
//
// ****************************************************************************

void
avtSimilarityTransformFilter::UpdateDataObjectInfo(void)
{
    //
    // The base class does some good work about setting extents.  Use that.
    //
    avtTransform::UpdateDataObjectInfo();

    if (atts.GetDoRotate())
    {
        double X = atts.GetRotateAxis()[0];
        double Y = atts.GetRotateAxis()[1];
        if (X != 0. || Y != 0.)
        {
            avtDataAttributes &inAtts = GetInput()->GetInfo().GetAttributes();
            avtDataAttributes &outAtts= GetOutput()->GetInfo().GetAttributes();
            if (inAtts.GetSpatialDimension() < 3)
            {
                //
                // Tell the new plot that it has spatial dimension 3.  Also
                // lie and say that its topological dimension is also 3.  This
                // will cause the facelist filter to execute.
                //
                outAtts.SetSpatialDimension(3);
                if (inAtts.GetTopologicalDimension() == 2)
                    outAtts.SetTopologicalDimension(3);

                double extents[6];
                vtkMatrix4x4 *t = GetTransform();

                //
                // Transform the extents into 3D.
                //
                if (inAtts.GetOriginalSpatialExtents()->HasExtents())
                {
                    avtExtents *newOriginalSpatial = new avtExtents(3);
                    inAtts.GetOriginalSpatialExtents()->CopyTo(extents);
                    extents[4] = extents[5] = 0.;
                    newOriginalSpatial->Set(extents);
                    *(outAtts.GetOriginalSpatialExtents()) = *newOriginalSpatial;
                    outAtts.GetOriginalSpatialExtents()->Transform(t);
                    delete newOriginalSpatial;
                }

                if (inAtts.GetThisProcsOriginalSpatialExtents()->HasExtents())
                {
                    avtExtents *newThisProcsOriginalSpatial = new avtExtents(3);
                    inAtts.GetThisProcsOriginalSpatialExtents()->CopyTo(extents);
                    extents[4] = extents[5] = 0.;
                    newThisProcsOriginalSpatial->Set(extents);
                    *(outAtts.GetThisProcsOriginalSpatialExtents()) = 
                                                     *newThisProcsOriginalSpatial;
                    outAtts.GetThisProcsOriginalSpatialExtents()->Transform(t);
                    delete newThisProcsOriginalSpatial;
                }

                if (inAtts.GetDesiredSpatialExtents()->HasExtents())
                {
                    avtExtents *newDesiredSpatial = new avtExtents(3);
                    inAtts.GetDesiredSpatialExtents()->CopyTo(extents);
                    extents[4] = extents[5] = 0.;
                    newDesiredSpatial->Set(extents);
                    *(outAtts.GetDesiredSpatialExtents()) = *newDesiredSpatial;
                    outAtts.GetDesiredSpatialExtents()->Transform(t);
                    delete newDesiredSpatial;
                }

                if (inAtts.GetActualSpatialExtents()->HasExtents())
                {
                    avtExtents *newActualSpatial = new avtExtents(3);
                    inAtts.GetActualSpatialExtents()->CopyTo(extents);
                    extents[4] = extents[5] = 0.;
                    newActualSpatial->Set(extents);
                    *(outAtts.GetActualSpatialExtents()) = *newActualSpatial;
                    outAtts.GetActualSpatialExtents()->Transform(t);
                    delete newActualSpatial;
                }
            }
        }
    }
}


// ****************************************************************************
//  Method: avtSimilarityTransformFilter::PostExecute
//
//  Purpose:
//      This is called to set the inverse transformation matrix in the output.  
//
//  Programmer: Kathleen Bonnell 
//  Creation:   April 10, 2003 
//
//  Modifications:
//    Kathleen Bonnell, Wed Jun  2 09:08:35 PDT 2004
//    Set both InvTransform and Transform in the output.
//
//    Kathleen Biagas, Mon Aug 15 14:17:39 PDT 2016
//    VTK-8, vtkMatrix4x4 no longer has operator[], use DeepCopy instead.
//
// ****************************************************************************

void
avtSimilarityTransformFilter::PostExecute()
{
    double m[16];
    vtkMatrix4x4::DeepCopy(m, invM);
    GetOutput()->GetInfo().GetAttributes().SetInvTransform(m);
    vtkMatrix4x4::DeepCopy(m, M);
    GetOutput()->GetInfo().GetAttributes().SetTransform(m);
}


// ****************************************************************************
//  Method: avtSimilarityTransformFilter::TransformData
//
//  Purpose:
//    Performs Translation and Scaling on the point-data scalars of a 
//    1D Rectilinear grid representing a curve.
//
//  Programmer: Kathleen Bonnell 
//  Creation:   March 28, 2008 
//
//  Modifications:
//    Kathleen Biagas, Wed Aug 22 16:00:03 MST 2012
//    Preserve data type of the scalars.
//
// ****************************************************************************

void
avtSimilarityTransformFilter::TransformData(vtkRectilinearGrid *rgrid)
{
    double yScaleVal = atts.GetScaleY();
    double yTransVal = atts.GetTranslateY();
    bool doScale = atts.GetDoScale() && yScaleVal != 1.;
    bool doTrans = atts.GetDoTranslate() && yTransVal != 0.;

    if (doScale || doTrans)
    {
        vtkDataArray *scalars = rgrid->GetPointData()->GetScalars();
        vtkDataArray *transScalars = scalars->NewInstance();
        transScalars->SetNumberOfTuples(scalars->GetNumberOfTuples());
        transScalars->SetName(scalars->GetName());
        transScalars->DeepCopy(scalars);
        for (vtkIdType i = 0; i < scalars->GetNumberOfTuples(); i++)
        {
            double val = transScalars->GetTuple1(i);
            if (doScale)
                val *= yScaleVal;
            if (doTrans)
                val += yTransVal;
            transScalars->SetTuple1(i, val);
        }
        rgrid->GetPointData()->SetScalars(transScalars);
        transScalars->Delete();
    }
}
